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ABSTRACT 

Two-dimensional (2D) hydrodynamical simulations of progenitor evolution of a 23 
Mq star, close to core collapse (in ~ 1 hour in ID), with simultaneously active C, Ne, 
O, and Si burning shells, are presented and contrasted to existing ID models (which 
are forced to be quasi-static). Pronounced asymmetries, and strong dynamical inter- 
actions between shells are seen in 2D. Although instigated by turbulence, the dynamic 
behavior proceeds to sufficiently large amplitudes that it couples to the nuclear burn- 
ing. Dramatic growth of low order modes is seen, as well as large deviations from 
spherical symmetry in the burning shells. The vigorous dynamics is more violent than 
that seen in e arlier burning stages in t he 3D simulations of a single cell in the oxygen 



burning shell (IMeakin fc Arnettll2007bl ). or in 2D simulations not including an active Si 



shell. Linear perturbative analysis does not captur e the chaotic behavior of turbulence 



(e.g., strange attractors such as that discovered by lLorena (|l963l )). and therefore badly 
underestimates the vigor of the instability. 

The limitations of ID and 2D models are discussed in detail. The 2D models, 
although flawed geometrically, represent a more realistic treatment of the relevant dy- 
namics than existing ID models, and present a dramatically different view of the stages 
of evolution prior to collapse. Implications for interpretation of SN1987A, abundances 
in young supernova remnants, pre-collapse outbursts, progenitor structure, neutron star 
kicks, and fallback are outlined. While 2D simulations provide new qualitative insight, 
fully 3D simulations are needed for a quantitative understanding of this stage of stellar 
evolution. The necessary properties of such simulations are delineated. 

Subject headings: stars: evolution, hydrodynamics: convection -turbulence, supernova: 
-remnants -core-collapse -nucleosynthesis 
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Introduction 



The first detailed calcula tions of the final, neutrino-cooled b urning stages prior to core collapse 
of massive sta rs were done b y lRakavv. Shaviv &; Zinamonl (|1967l ). with simplified energy generation 
rates, and by lArnettl ()1968l ) using a small (24 species) network. Because of the extreme demands 
upon computer resources then available, simplified energy generation rates were used in subsequent 
calculations; see lArnettl (119961 ) for references to the early work. This work showed two issues which 
have not yet be resolved: (1) the definition of convective zo ne boundaries i nvolves incomplete physics 
(by construction mixing length theory ignores gradients, Spiegell (1971)), and (2) the stability of 
the burning in a convective region has not been demonstrated, only assumed. While thermal 
instability (thermal stability implies a global balance betwe en nuclear heating and neutrino cooling 
in the convective zone) was considered (|Arnettlll972al . ll996l )). dynamical instability (i.e., instability 
related to fluid flow) was not. In this paper we will show that both these issues are significant, and 
that they can be resolved by 3D numerical simulations and theory. 

Almost all previous progenitor models for core collapse have focused on thermal behavior and 
quasi-static mixing, which are described by the evolution of the temperature and the composition 
variables. The dynamic behavior of the stellar plasma includes and is dominated by the veloc- 
ity fields, which not only imply mixing, but also possible change in the stellar structure. The 
star is not necessarily a quasi-static object, but may have significant fluctuations in its variables. 
The dynamical behavior, found here for simultaneous, multi-shell burning, will drive entrainment 
(jMeakin fc Arnettl l2007bl ) at convective shell boundaries, changing the nucleosynthesis yields and 
the size of the Fe core at collapse. 

The presupernova structure of a massive star consists of a core, mantle, and envelope. The 
envelope is extended, composed of H and He, and may have been removed prior to core collapse 
by wind-driven mass loss, or tidal stripping by a companion. The mantle is composed of burning 
shells of C, Ne, O and Si; these shells are convective, interact nonlinearly, contain most of the 
nucleosynthesis products ejected, and may smother a neutrino-driven explosion. The core is com- 
posed of Fe-peak nuclei, and its mass is determined by its entropy and electron fraction. Lower 
entropy and electron fraction give smaller cores, which are easier to explode by neutrino transport 
mechanisms. Core collapse mechanisms for explosion are sensitive to core mass, mantle density, and 
rotation; all are sensitive to the treatment of turbulence. The simulations described here involve 
the simultaneous action of C, Ne, O and Si burning shells. The oxygen shell is special, because 
(1) formation of electron-positron pairs softens the equation of state, aiding the formation of large 
amplitude motion, (2) the large abundance of oxygen and its relatively large energy release per 
unit mass provide a large thermonuclear energy reservoir, and (3) oxygen burning, unlike silicon 
burning, has little negative feedback from quasi-equilibrium to damp flashes (see below). 

In Section 2 we discuss the historical context of progenitor models of core collapse supernovae, 
focusing on issues of mixing, causes of time-dependence and multi-dimensionality, prospects for 
development of better computational tools, and the differences between 2D and 3D simulations. 
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In Section 3 we describe our 2D simulations of progenitor evolution with multiple, simultaneously- 
active, burning shells (C, Ne, 0, Si), and discuss some new phenomena which appear. In Sec- 
tion 4 we summarize the implications for evolution prior to core collapse. In Section 5 we focus 
on several problems in observational astrophysics which need to be reconsidered in light of our 
results. In Section 6 we summarize the major conclusions, and outline the necessary features of 
new 3D simulations which will be required to quantitatively resolve the issues presented by the 2D 
simulations. 



2. Brief Historical review of ID, 2D, and 3D Models 

2.1. ID Models: Mixing 

By the term advection we mean the transport of a parcel of matter by a large-scale flow; by 
diffusion we mean the transport of a parcel of matter by a random walk of small-scale motions 
(often microscopic motion of ions in the stellar plasma) . The mass (baryon) flux due to advection 
is 

F m = up, (1) 

where u is the fluid velocity vector and p is the mass (baryon) density. The flux of any scalar 
variable is related to this flux by a factor of the density of the variable per unit mass (baryon). For 
example, if the mass (baryon) fraction of nuclear species i is Xj, the flux of this species is 

Fi = u P X t = F m X t . (2) 

The rate of change due to such fluxes involves a divergence, as in the continuity equation, 

dp/dt = -V • F m . (3) 

Thus, ID advection involves a first order spatial differential operator. 

If the considered volume contains heterogeneous (turbulent) matter, the flux may vary over 
the corresponding surface. If that volume is a zone in a stellar evolution computation, then the 
fluxes must be averaged over the heterogeneity, and some knowledge (or assumption) regarding the 
smaller scale structure is required. Consider the simple example of only an inflow and an outflow, 
and variation only in the vertical direction. The net flux for species Xi is then 

FnetiP^i) — {FmXi)out {.-^m-^-i)in- (4) 

Suppose that F out = Fi n = pu, so F out — Fi n = 0, which automatically satisfies mass (baryon) 
conservation, so that 

F net (Xi) = -pu {Xi) out - (Xi) in . (5) 
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In this case the flux (of composition i) is proportional to the negative of the difference in abundance 
in the up and down flowa[j]. 



For diffusion, the number flux follows from Fick's law, 

1 



(6) 



and is proportional to the gradient of a number density iVj. Here A is a mean- free-path and v is 
a speed of diffusing particles. The number density for species i may be written as iVj = pAXi/Ai, 
where A = l/m amu is Avogadro's number (the inverse of the atomic mass unit) and Ai is the 
number of nucleons in species i. The mass (baryon) flux is FdpQ) = Aim amu Fd(Xi). Thus, 



F d (X, 



l -\vV{pXi). 



(7) 



If we consider the divergence of the flux, diffusion implies a second order spatial differential operator, 
in contrast to advection which, as we saw, implies first order. This is a fundamental mathematical 
difference. 



Weaver. Zimmerman fc Wooslevl (|1978l ) ; IWooslev. Heger. Weaver! (|2002l ) introduce an effec- 
tive diffusion coefficient D to model convection: 



dt 



d_ 

dm 



(Airr 2 fp 2 D 



dm 



(8) 



In the case of a region unstable to convection by the L edoux criterion, they take D a = ty.l /3, where 
v c is the velocity estimated by mixing- length theory (jBohm-Vitensd Il958l ; IClavtonl Il983l ) and i is 
the mixing length. Thus if 



Xv = iv c 



(9) 



then Equation [7J is recovered. The use of Equation [7J requires that A <C Ar, where Ar is a zone 
size. For a turbulent cascade, I ^> Ar, which is inconsistent with the previous requirement. The 
treatment of convection as diffusion is essentially an algorithmic interpolation procedure, but has an 
intrinsic contradiction in physics. This arises because turbulent flow has two facets: a flow at large 
scale I which does most of the transport, and a flow at small scale A which does the dissipation, and 
£ S> A. The diffusion approximation might be used for the small scale flow, but that is irrelevant 
to th e transport problem, which is dominated by the large scale flow (jArnett. Meakin. &; Young 
2009h . 



Not ice the rough similarity between Equation [5l which represents the underlying fluid dy- 
namics ([Landau &: Lifshitzlll959l ). and Equation [7J which represents the proposed approximation. 
Taking the density outside the differential operator and using a finite difference representation of 
Equation [7J 



Fd(Xi) 



—Xvp/Ar 



{X, 



ilout 



(Xi 



(10) 



A more realistic and relevant case for stellar turbulent convecction would allow different speeds in up and down 
flows; that complication is not necessary here. 
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For this to be consistent with ID advection (Equation [5]) we must have a diffusion rate that is depen- 
dent upon the z oning! While perhaps usef ul algorithmically, this is not clear conceptually; see also 
the discussion of Maeder Sz Mevnet ( 200ol ). who also express doubts concerning the approximation 
of advection by diffusion. 

Mathematically, diffusion is a procedure which maximally smooths gradients, so that a more 
realistic procedure may be expected to exhibit less smoothing. While based on poor physics, the 
diffusion model of convection had the virtue that it allowed numerical prediction of yields. 

M uch effort has recently gone into extending stella r evolution to include ro t ational effects 



( Zahn 



1992 



Chabov er. Demarque. fc Pinsonneaultl Il995l : iMaeder Meynetl I200Q : iTassoull l2000l : 



Woosley. Heger. Weaverl 120021 ) . Evolution of a rotating star will develop differential rotation in 
general, and shear flow. Similarly, deceleration of convective plumes will develop shear flow. In stars 
both flows will be turbulent. Convection and rotation have an unde rlying similarity not reflected 
in stellar evolution theory. The review of IMaeder &: Mevnetl (|2000l ) gives a clear presentation of 
the various approximations involved in reducing the full fluid dynamic equations to a simpler, 
more easily solvable se t. It is now possible to test these approximation s by numerical simulation 
i n both shearing box (lArnett Meakin 2009j; IStone fe: Gardinerl 120101 ) and whol e star domain s 
dBrun fc Palacioslhood : brown, et al.l hoirtn" see also the theoretical developments of lBalbusI ([2009); 
Balbus fc Weiss] (|2010l ). A theoretically sound approach must treat the shear from differential 
rotation and the s hear from con vective plumes on a consistent basis, reflecting their underlying 
physical similarity (|Turnerlll973l ). 



2.2. Time Variation 



The e-mechanism. The e-mechanism ( 



edouxj Il94ll . Il958l ) for driving stellar pulsations by nu- 



clear burning was examined by lArnettl (119771) for Si burning in I D geo metry, using the simple 



energy generation rate proposed by iBodanskv. Clayton h Fowlerl (j 19681 ). In th is case, explic it 
ID hydrodynamic simulations gave nuclear-energized pulsations (see Figure 6 in lArnettl (|1977l )). 
There were two high frequency modes (period 0.1 seconds and 1 second) and a lower frequency 
convective mode (turnover time ~ 20 seconds). The intermediate frequency mode was related to 
the Si flash. Computation with a realistic nuclear network subsequently showed that the highest 
frequency ("acoustic") pulsations would be strongly damped due to the quasi-equilibrium nature 
of Si burning; an increase in temperature gave an increase in free alpha particles (and neutrons 
and protons), which required energy, and resisted the increase in temperature (with an almost 180 
degree phase shift, meaning that there is strong negative feedback to resist changes). However, this 
damping process does not apply to O, Ne or C burning, which have little or no quasi-equilibrium 
behavior, and in principle could drive pulsations more vigorously. 



The_r- mechanism. lArnett k, Meakinl (j2011bl ). in study of the 3D simulations of lMeakin Sz Arnett 
(|2007bl ) , found that the bursts in turbulent kinetic energy were related to similar fluctuations in the 
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Lor end (|1963l ) model of a convective roll; this is the now-famous strange attractor. This instability 
mechanism, called the r-mechanism for "turbulence" , is expected to be a general property of stellar 
convection zones, and probably is the cause of th e fluctu ations in luminosity observed in irregular 
variables (e.g., Betelgeuse, see lArnett & Meakinl (|2011al lbh. This process is independent of the 
temperature and density dependence of the thermonuclear heating, and thus is distinctly different 
from the e-mechanism. Unlike the e-mechanism, the r-mechanism is not a linear instability, but is 
inherently nonlinear. 

In appropriate circumstances however, the r-mechanism may couple to nuclear burning, mak- 
ing pulsations more violent, giving a more complex, combined e+r mecha nism. The r mechanism in- 
volves a non-linear instability, unlike the linear instabilities disc ussed bv lGoldreich. Lai fe: Sahrling 
(| 19971 ). In a detailed analysis. iMurphy. Burrows fe Hegerl (120041 ) solve the linear perturbation equa- 
tions for several massive stars prior to core collapse, and based upon the e mechanism alone, find 
inadequate driving to cause such violent behavior as is actually shown in the non-linear simulations 
presented below. 



The analysis of lMurphv. Burrows Hegerl (|2004l ) is perfectly correct within framework of the 
linear assumption, but the assumption fails. The linear appr oximati o n is n ot valid in this case 
because turbulence is not a linear perturbation of the system. Lorenz (jl963l ) showed that chaotic 
behavior in a convective roll is due to the nonlinear interaction between temperature gradients 
(both vertical and horizontal) and convective speed. Solutions which are initially close to each 
other will diverge exponentially as time passes. Thus the interplay between pulsation and turbulent 
convection is not captured by traditional linear perturbation analysis of stellar pulsations. T his is 
an in t eresting theoretica l result, suggesting that linear perturbative methods for pulsations (|Cox 



1980 



Unno. et al.lll989l ) may require reexamination when turbulent convection is important (as 



those authors feared). 



2.3. History of Multi-dimensional Progenitor Models 



There have been relatively few multi-dimensional simulations of core colla pse progenitors, but 



there has been a rich context of efforts on turbulence and stellar hydrodynamics 



2006; Mocak. et al. 



ulations by 



2008 



Brownind (|200 



Porter k, Woodward 



1994 l2Q00l:rFrvxell. Miiller fc Arnettlll989h . on the helium flash (lDeupreelll984l . ll996l . l200.ilDearborn. LattanzioTfc 



, \j X II -1 7 1 i ij I 

2 0091 . boiol ). and on turbulent MHD with rotation (e.g., m-dwarf sim- 
as well as extensive work on the Sun with the ASH code ( e.g., 



Brown, et al"1(|201oh ) and on stellar atmospheres pioneered by Nordlund and Stein (e.g.. lNordlund. Stein, k, Asplun 
(|2009l )). This context has speeded the develo pment of too ls and helped determine their reliability. 
The first effort at simulating oxygen burning (jArnettll 19941 ) helped define the shell burning problem 
with regard to required resolution, but suffered from the use of sectors so narrow in angle that 
boundary effects affected the flow. 



Table [T] summarizes aspects of early 2D simulations of progenitor models in comparison to 
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the present work. The first 2D hydrodynamic simulations of core-collapse progenitors (oxygen 
burning) showed striking new phenomena: mixing beyond formally stabl e boundaries, hot spot s 
of burning du e to C 12 entrainment, a nd inhomogeneity in neutron excess ( Bazan fc ArnettJ 11994 ). 
Further work (jBazan fc Arnettl Il998l ) confirmed that convection was too dynamic to be well rep- 
resented by diffusion- like algorithms, that large density perturbations (8%) formed at convective 
bou ndaries, and that gravity waves were vigorously generated by the flow. Extension to Si burn- 
ing (jBazan &i Arnettl 1 1 99 7bl ) with a 123 isotope network showed similar highly dynamic behavior 
and significant inhom ogeneity in neutron excess. With an entirely different 2D code (VU LCAN), 



Asida Arnettl (|2000l ) extended the evolution of the oxygen shell of lBazan Arnettl (|l998l ) to later 



times, and discovered that the extensive wave generation at the convective boundaries induced a 
slow mixing in the bounding non-convective regions. 

Kuhlen. Woosley. Glatzmaierl (hood) investigated shell oxygen burning in 3D with an anelas- 
tic hydrodynamics code (jGlatzmaierl 11984 ) , and found small density and pressure perturbations 
(less than 1%). The boundary conditions were impermeable and stress- free, so that convective 
overshoot could not be studied. They concluded that, contrary to previous work listed above (done 
with 2D compressible codes), the behavior was not ver y dynamic and could be described by t he 
MLT algorithms used in ID stellar evolution codes (e.g.. IWeaver. Zimmerman Sz Wooslevl (|1978l )). 



Meakin &i Arnettl (|2007al ). using 3D compressible hydrody namics, showed that the differenc e s were 



due to the choice of (unrealistic) boundary conditions that iKuhlen. Wooslev. fc Glatzmaierl (|2003l ) 
used. Inside the convection zone, away from the boundaries, the Gl atzmaier code gave res ults in 
good agreement with the compressible code. However, as stressed in lBazan &: Arnettl ()1998l ). fluid 
boundaries al low surface waves to build to lar g e amp litudes (8p/p ~ 0.1), so that the hard bound- 
aries used in IKuhlen. Woosley. Glatzmaierl (120031 ) distorted the physics. The good agreement 
between the anelastic and the compressible solution within the convection zone, and the agreement 
between the stable layer dynamics given by the compressible fluid code and analytic solutions to 
the non-radial wave equation, indicate that the compressible hydrodynamic techniques are robust 
for this problem, even for Mach numbers below 0.01. An anelastic code with the correct boundary 
conditions should give the same result; the flow is subsonic. 

With the PROMPI code (|Meaki nlhoOfil ) . several new series of computations were done. iMeakin &: Arnettl 

(|2007al lbh presented the first 3D calculations of the phase of shell O burning with fu ll physics (i.e., 



Meakin & Arnett 



compr essibility, nuclear network, real equation of state, appropriate boundaries, etc., 

calculated in 2D the o xygen burning shell, and both C and Ne burning shells above this. 
Here we present simulations (jMeakinl 120061 ) with similar microphysics, extended to include the 
silicon burning shell as well, so that C, Ne, O and Si burning occur on the grid, but only in 2D. 



2.4. Future Prospects: Beyond 2D 



Progress will not be a simple progression reflect i ng the growth of computational resources, but 
also of theoretical understanding. lArnett &z Meakinl (|2011bl ) have shown a connection between the 
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Loren z model of a conve c tive ro ll and the bursty pulses in turbulent kinetic energy seen originally 
in the lMeakin &: Arnettl (|2007bl ) simulations, which seem to have a similar strange attractor. The 
Lorenz model is a low order (3 variable) dynamical system, and its physical identification with the 
3D simulations suggests how we may project the essence of the 3D solutions onto ID for stellar 
evolution and dynamics ("321D"). This will give an immediate improvement over MLT as well 
as better initial models for numerical simulation; see below. The physical basis is strengthened 
mathematically by use of t h e Ka rhunen-Loeve or proper orthogonal decomposition (POD; se e 
Holmes. Lumlev. &: Berkoozl (|1996l )) of the 3D numerical data set of iMeakin &: Arnettl (|2007bl ) . 
Preliminary results show that roughly half of the turbulent kinetic energy is in the single lowest 
order empirical eigenmode, supporting the idea that low order dynamical systems may be used to 
describe the complexities of time dependent turbulent flow in stellar convection zones. 

The dynamical phases prior to core collapse may not attain a statistically averaged state, so 
that these methods (321D and KL decomposition), while promising for earlier evolu tion, may no t 
be the optimum tools for the pre-collapse itself, where large eddy simulations (ILES. iBorid (|2007l )) 
in 3D are needed. However, use of the theoretical tools (321D and KL decomposition) provides a 
means of estimating the range of fluctuations about an given ILES simulation, which may be tested 
insofar as computational resources allow, for other initial conditions. 



2.5. Comparison of 2D and 3D Simulations 



Meakin Arnettl (|2007bl ) show a precise comparison of two simulations, which use the same 
microphysics, initial model and code, and differ only in dimensionality, which changes from 2D to 
3D. The simulations are discussed there as "core convection" (see Figure 4 therein). The topology 
of the convective flow is significantly different between 2D and 3D models: the 3D convective flow 
is dominated by small plumes and eddies while the 2D flow is much more laminar, and dominated 
by large vortices ("cyclones") which span the depth of the layer. The 2D vortices trap material 
which is slow to mix with surrounding matter; in 3D the flows become unstable and matter mixes 
more completely. The wave motions in the stable layers do not have an identical morphology in 
2D and 3D, and the velocity amplitudes are much larger in 2D. 

The differing behavior is due to the constraint of geometry, and the law of angular momentum 
conservation, which forces the formation of large cyclonic patterns in 2D that are unstable in 3D. 
The turbulent cascade moves from small scale to large (cyclones) in 2D, but from large to small 
(Kolmogorov cascade) in 3D. Cyclonic behavior at large scales is physically reasonable for the 
Earth's atmosphere, for example, which because of its restriction in height, is approximately two- 
dimensiona@. In stars there is no physical constraint to enforce 2D motion, so that 2D simulations 



2 The density scale height in the vertical direction is small compared to the width of a typical cyclonic system seen 
on the evening news. Oxygen would be required in an unpressurized aircraft at 35,000 feet, which is a measure of the 
"height" of the atmosphere. This is less that one percent of the width of a large cyclonic system. This is a very flat 



- 9 - 



of stars are not realistic, merely computationally less expensive. 



In lMeakin fc Arnett] (|2007bl ) , 3D simulations were done for a single cell in the O shell; the com- 
putational demands for simultaneous multi-shell burning are more extreme (but almost feasible). 
As a first step we present 2D simulations, which though flawed, are instructive. We describe t he 
first extended results for four simultaneous burning shells (Si burning included, Meakin (j2006l )). 
and interpret them with the benefit of extensive analysis of 3D O burning results. Taken with 
Figure 4 from iMeakin &; Arnettl (|2007bl ). they provide clues about the nature of the true behavior 
to be expected in 3D. 



3. 2D Simulations of Multiple Simultaneously Burning Shells 

3.1. Starting Point 

The initial conditions are a fundamental problem for numerical simulations which explore the 
development of instabilities to large amplitude. Subtle biases in the initial state might make the 
subsequent development misleading. Use of a ID initial model, with no reliable description of the 
turbulent velocity field or the extent and position of the boundaries of the convection, is a cause 
for worry, but is the best we can do at present. 

The 2D si mulations started from a mo del obtained from a sequence by P. A. Young (private 
communication; iFryer. Young, et al. I (|2008l ) have evolved this model to core collapse in ID, and to 
3D explosions). The initial state was a ID model of a 23 Mq star, mapped to 2D; it is at the latest 
stages of evolution, about one hour prior to core collapse. Turbulent convection developed from very 
small perturbations (~ 10~ 3 or less, due to mapping onto a different grid) in the unstable regions. 
The convection rapidly evolved to a dynamic state with far larger fluctuations , indepe ndent of the 
initial perturbations. The model was similar to that used in lMeakin &: Arnettl (|2007bl ) . except the 
computational volume is deeper, reaching down past the Si burning shell, the aspect ratio is larger 
(a full quadrant was calculated), the abundance gradients were smoother (more like a diffusion 
approximation for convection; see § 2.1), and the onset of core collapse was much nearer (~ 1 
hour). 

The computational domain had an inner boundary representing the Fe core and extended 
well beyond the active burning regions to the edge of t he He core. The bounda r y cond itions in 

The 



2007. 



3Q. 



angle were periodic and in r adius were reflecting, as in IMeakin Arnettl ((2006, 
equation of state was that of iTimmes &: Swestvl (|2000l ). which accurately describes the effects of 
partial relati yistic degeneracy of elec trons and of the formation of electron-positron pairs, and is 
very similar (jTimmes &: Arnettl I1999 ) to the equation of state used in previous simulations listed 
above. 



domain, and with rotation, is dominated by geostrophic (i.e., 2D) flow patterns at the large scales. 
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Table 1. Two-dimensional Progenitor Simulations 



Reference 


a 


b 


c 


d 


e 


zoning 


256x64 


256x64 


172x60 


800x320 


800x320 


code 


PROMETHEUS 


PROMETHEUS 


VULCAN 


PROMPI 


PROMPI 


eos f 


wda 


wda 


wda 


TS 


TS 


network 


12 


123 


12 


25 


37 


burning 





Si 





C,Ne,0 


C,Ne,0,Si 


core 


Si 


Si 


Si 


Si 


Fe 


duration(s) 


300 


200 


1200 


2500 


600 



Bazan fe Arnettl (1998, 1994) 



Bazan fe Arnettl (|l997bl b a small inert spherical boundary surrounded by Si. 



Asida fe Arnettl (|2000h 



MMeakin fe Arnettl (|2006T ). energy release verified against a 177 nuclei network 
Meakin ( 2006h and this paper, energy release verified against a 177 nuclei network 



See lTimmes fe Arnettl (|l999 ) and lTimmes fe Swestvl (|2000h 
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Fig. 1. — Thermonuclear reaction network (37 nuclei) used for burning in C, Ne, O, and Si shells. 
Each box represents a nucleus; see text for details. 



- 11 - 



Figure [T] shows the nuclear reaction network used, which contained 38 species (37 isotopes 
plus electrons). To deal with increasing neutron excess it was extended to Ni 62 , which corresponds 
to Z/A ~ 0.45. During Si burning, the nuclei having Z > 22 approach a local nuclear statistical 
equilibrium (they become a quasi-equilibrium group). The results of this network were compared 
to those of a 177 isotope network; it reproduced both the energy generation and the increase in 
neutron excess predicted by the larger n etwork ( see di scussion of nucleosynthesis and of increase in 
neutron excess during silicon burning in I Arnettl (1996)). The nuclear burning was directly coupled 
to the fluid flow by the method of operator splitting. 



3.2. Results 

Figure [2] shows the structure of a quadrant of the core of the 23 M star, with an iron core 
(shown as a white semi-circle) in the center. The computational domain includes burning shells of 
Si, O, Ne and C, in order of increasing radius. The left quadrant displays abundance contours of 
Si 28 (dark blue is high abundance, white is low). The inner light region is a burning shell which is 
depleting its Si. Above this is a dark ring which is unburned Si. Further out is a medium blue layer 
which contains the O burning shell, in which Si is being produced. At the top of this layer are seen 
streams of very light blue, denoting entrainment of matter which has not been contaminated by 
oxygen burning (mostly C and Ne). Finally there is a very light blue layer from which the streams 
came, and which has no enhancement of Si above its original value. 

The right quadrant shows contours of energy generation rate, in units of 10 13 erg/s. Note 
that the second quadrant is presented as a rotation about the vertical axis; this helps identify 
corresponding matter in the two different variables which are mirror images. The inner ring is 
again the Si burning shell. It has both strong heating (yellow) and strong cooling (blue) at the 
same radius, that is, the burning does not possess spherical symmetry. This is enclosed by a green 
ring, the Si rich layer, which has milder neutrino cooling, and is no longer spherically symmetric 
because of "hot spots" of burning in entrained (descending) plumes. Beyond this is a ring of red 
and yellow which is highly dynamic: the O burning shell itself. Finally, above this wisps and plumes 
of heating are beginning to be seen; these are due to C and Ne burning in entrained matter which 
is rich in these fuels. 



These models have a low Damkohler number (jDamkohlerl ll94G ) . D a = Turbulence/ Tread ^ 1) 
where the time scale for turbulent flow is Turbulence an d the reaction time is T reac t. In this limit 
a complex mixing model is not needed to describe the burning, unlike burning fronts in type la 
supernova models which require a more complicated description. Here fuel is transported into 
regions of higher pressure, compresses and heats, burns, expands and cools, and is buoyantly 
transported back to lower pressure. Changes in composition due to turbulent mixing are much 
faster than those due to nuclear burning, so that a sub-grid flame model is not required. 



The initial model was strictly spherically symmetric, and had to develop a self-consistent 
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Fig. 2. — A snapshot of the structure of C, Ne, O, and Si shells surrounding the Fe-core of a pre- 
collapse progenitor of 23 M star, about 500 seconds after the constraint of spherical symmetry 
has been removed. The left side shows the abundance of Si 28 while the right side shows the net 
energy generation rate. 
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Fig. 3. — Snapshots of the structure of C, Ne, O, and Si shells surrounding the Fe-core of a pre- 
collapse progenitor of 23 Mq star. Three different times are shown, tt = 0, 61, 83 seconds (from 
top, s, to bottom, 83 s) after our fiducial model (see text). The left panels (blue) show abundance 
of Si 28 , while the right panels show energy generation rate and convective speed, respectively. 
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turbulent velocity field to carry the heat from nuclear burning away from the burning regions. 
During this mild transient phase, asymmetries begin to develop slowly. Figure [2] shows the structure 
after this development, and at the beginning of significant deviation from the usually assumed 
spherical symmetry. Because of the time needed for the initial spherical model to develop a realistic 
and consistent convective flux, there is ambiguity regarding the "initial time" ; we simply define a 
"fiducial time" (tf = t — to where to ~ 345 seconds into the 2D simulations) at which the model is 
still fairly spherical but has a realistic convective flux, and we count elapsed time after that. 

Figure [3] shows three snapshots of the structure at different times after the fiducial time: tf = 
0, 61, and 83 seconds. The left column represents the same variables as shown in Figure 1. The 
right column shows the Si 28 abundance as before in the left panel, but in the right panel the energy 
generation rate is replaced by the turbulent speed in units of 10 7 cm/s. 

Distortions in the O burning shell are obvious. The equation of state in this region is affected 
strongly by the thermal production of an equilibrium abundance of electron-positron pairs, so that 
the effective adiabatic exponent drops below{§ Ti ~ 4/3. Similarly T4 = 1 + E/PV ~ 4/3; This 
means the local contribution to the gravitational binding energ y, which is proporti onal to T4 — 4/3, 



is small. This is a common property of oxygen burning in stars lArnettl (|1968l . Il996l ). The restoring 



force for stable stratification is weak, allowing large amplitude distortions with little energy cost. 

The decrease in T\ and T4 are due to the need to provide the rest mass for newly formed 
electron-positron pairs. At low temperature, kT <C 2m e c 2 , the number density of pairs relative to 
the charge density of ions decreases exponentially with decreasing temperature, and is negligible. 
As temperature approaches T ~ 2 x 10 9 degrees Kelvin, the effect on the equation of state becomes 
largest. At higher temperatures, the increase in mass of pairs is less relative to the thermal energy 
kT, so that these gammas approach that of an extreme relativistic gas, Fi = T4 = 4/3. Oxygen 
burning (O 16 fusion) in stars occurs at T ~ 2 x 10 9 de grees Kelvin, so that this burning stage is 
most influenced by the effects of electron-positron pair production on the equation of state. 

Consider first the left column in Figure [3j The top panel (t=0) is relatively symmetric, but 
as time passes, the middle and lower panel show increasing distortion, especially visible at the 
interface between the outer, light blue layer and the middle, medium blue layer. The streams of 
light blue inside the medium blue represents entrainment of matter with little Si, that is, C and 
Ne fuel. This corresponds to the flame structures seen in the right side of the left column, which 
are due to C and Ne burning (note similarity in shapes on left and right sides in the left column). 
Similar entrainment is occurring at the top of the Si burning convective region; the outer edge of 
the light blue inner ring is rippled due to bursts of burning. The amplitude of these variations is 
smaller due to the stiffer equation of state here. 

The right side of the right column shows the turbulent convective speed. The large structures 



See Table 5 in lArnett. Meakin. fc Younj (120091 ). and recall that X7 ad = (T 2 - 1)/T 2 ~ 1/4 across the whole O 



burning convective zone. 
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are the oxygen burning convective zone. A smaller convection zone may be seen surrounding the 
Fe core, due to the Si burning shell. The C-Ne layer, lying outside the oxygen burning shell, 
illustrated the effect of a low-order mode. In the top and bottom panels, there is little motion, 
while in the middle panel the amplitude of the motion is near maximum. The lighter red areas at 
about 30 degrees and 70 degrees from vertical correspond to nodes in the modal velocity. Because 
of symmetry about the equator (horizontal) we have four nodes in 180 degrees, or an I = 4 mode 
being dominant. Odd values of t are suppressed by the domain size and symmetry imposed by our 
boundary condition, so this is the lowest order possible in this simulation; it has a period of about 
60 seconds, but is mixed with other, weaker modes. A movie of the simulation shows a dramatic 
change as this mode turns the speed on and off as we move from top to middle to bottom. 

Figure H] shows the same variables at t = 115, 247, and 307 seconds, after several, increasingly 
vigorous "sloshes". The distortion in the O shell has grown, with large amplitude motions whose 
wave- form is characteristic of g-modes. Entrainment of C and Ne increases. Fluctuations in Si 
burning become more vigorous and distort the Si layer although quasi-equilibrium damps explosive 
excursions. The Si shell burning occurs in dynamically form i ng and disrupting cells, qualitatively 



similar in nature to the 3D cell studied in lArnett fe Meakinl (|2011bl ) . A pink tinge indicates weak 



but widespread burning of carbon in the outer layers of the computational domain. 

At tf = 307 seconds what is best described as an "eruption" is occurring. The distortion in 
the O shell (medium blue bulge) continues to grow; it is due to strong wave motion, powered by 
the nuclear burning. In the bottom panel, the thickness of the O shell convective layer varies by 
more than a factor of three as a function of angle; it is thin at the equator and thick at a 45 degree 
angle. This corresponds to an eruption at that angle. Entrainment is correspondingly increased, 
with consequent burning. Waves "slosh" back and forth. The computational domain and boundary 
conditions restrict the lowest order modes. The Si burning shell is affected by the behavior of the 
O shell. There appears to be no evidence of a slacking in growth of the dynamic behavior, and 
extensive mixing is occurring. We ended the simulation because the extent of mixing is so great that 
it reaches the grid boundaries, and the simulation domain becomes inadequate. The assumption 
of spherically symmetric structure has clearly failed, as has the assumption of quasi-hydrostatic 
structure. 



4. Implications 

Figures 3-4 show the breakdown of the assumption of spherical symmetry, which has been 
the basis of supernova progenitor models, and much of the interpretation of observational data of 
supernovae and young supernova remnants. The t = 4 spherical harmonic is dominant because the 
computational domain had only one quadrant and periodic boundaries; a full hemisphere would 
have allowed the I = 1 mode to appear. 

These simulations, of a stage only one hour before core collapse, show more violent behavior 
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Fig. 4. — Structure of C, Ne, O, and Si shells surrounding the Fe-core of a pre-collapse progenitor 
of 23 Mq star. Snapshots at times tt = 115, 247, and 307 seconds (top to bottom) after our 
fiducial model (see text). The format is the same as in Figure O The eruption has become strongly 
non-linear, as the bottom panels show. 
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than previous 2D simulations of multiple (CNeO) shells (jMeakin fe Arnettll2006l ) , or 3D simulations 
with a single burning shell (jMeakin fe Arnettl l2007bl ) . They become more violent in less time; 
compare durations in Table 1. There are several reasons for this behavior. (1) This stage is late, 
only an hour prior to collapse. (2) The 3D O shell simulations had a smaller computational domain, 
so that low order modes were restricted by periodic boundary conditions. (3) The difference between 
the 2D three shell (CNeO) and four shell (CNeOSi) simulations seems to be due to an additional 
interaction between the O and the Si shells, which are both active. 

Entrainment of new fuel causes mixing, which is countered by heating from burning, which 
changes entropy deficits in down-flows to entropy excesses, and halts the down-flow motion. This 
gives a natural layering mechanism for the highly dynamic system, a dynamic layering. Obviously 
the compositional structure, and the predicted yields, depend upon this effect, which has been 
ignored (or treated in a diffusion rather than advection algorithm) in all ID models of progenitor 
evolution. Inter-zonal mixing also would have an impact on yields. 

Strong wave genera tion is observed. Such waves may become compressional (mixed mode, 
Meakin fe Arnettl (|2006l )) as they propagate into the strong density gradient. The waves will 
dissipate in non-convective regions, causing heating and slow mixing there, and to the extent that 
the wave heating is faster than radiative diffusion (which is very slow), expansion will occur. These 
effects are large enough to be seen in the simulations, but need further quantification to determine 
their relative importance. 

The Fe core (here a static boundary condition) will contract, giving increasingly dynamic 
behavior. Such a dynamic approach to core-collapse has not been investigated in multi-dimensional, 
multi-shell simulations, but may have interesting consequences (see below). The duration of the 
2D CNeOSi simulation was ~ 600 seconds, whil e the time fo r core collapse, neutrino trapping, 
and rebound shock are together a few seconds (I Arnettl 1 19961 ). The observed diffusion time for 
neutrinos from SN1987A was comparable, also a few seconds. Because time for core collapse is fast 
in comparison to the period associated with O shell dynamics, the shell structure may be caught 
in a distorted state by the explosion shock, giving a no n-spherical remn a nt ev e n if t he explosion 
shock were perfectly spherical (which it is unlikely to be: iKifonidis. et al. I (|2003l . 120061 )). 



The 2D simul ations show much more active dy namical behavior than suggested by linear per- 
turbation analysis (IMurphy. Burrows fe Hegerll2004r) . In ^2.2l this was traced back to the treatment 



showed that convection has an intrinsic 



of turbulent convection in the linear analysis, 
nonlinear instability; it arises from advection, which gives product terms (XY and XZ in his no- 
tation); these are produc ts of the velocity amplitu de with the horizontal and vertical temperature 
fluctuations respectively ( Arnett fe Meakin 201 lbn. This provides an ex plicit example in which lin- 



ear perturbative methods for pulsations (jCox 



1980 



Unno. et al.l Il989l ) require modification when 



turbulent convection is present. Study of low order dynamical models, such as lLorenzl ()1963l ). will 
provide insight into the nature of this general problem. 



It is unlikely that ID progenitor models are realistic; in addition to unlikely geometry, they 
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ignore the vigorous dynamic behavior, which becomes manifest in 2D and 3D simulations, in which 
it is not forbidden as it is in ID. There may be eruptions prior to core collapse, there will be large 
amplitude distortions away from spherical, "onion" -skin structure, and there may be modifications 
of the supernova shock by explosive oxygen burning. Not only can explosion shocks be non-spherical, 
but the progenitor mass they propagate through can be asymmetric as well. The distortions in the 
progenitor and those induced by the explosion shock may leave an imprint on the abundances in 
the supernova and its ejecta. Of course, rotation will have its own effects which we will have to 
disentangle. 



Some Issues to Reconsider 



Progenitor Structure and Fall-back. Mechanisms of explosion and fall-back are all predi- 
cated upon the validity of the detaile d structure o f the ID progenitor models. The old problem 
of non-explosion of supernova mo dels (IArnettlll996l) has been alleviated by multi-d i mensional sim- 



ulatio n s of collapsing cores (e.g., iBurrows. Livne. Dessart. Ott. &: Murphvl (|2006l ); iBruen. et al. 



(|2009l ): IWongwathanarat. Janka. &: Mulled (j2010l )). More realistic initial models will bring further 
changes. For example, the explosion and remnant formation could be modified if an eruption oc- 
curred prior to and during core collapse. Similarly, expansion of the mantle surrounding the Fe core 
could be caused by dynamic burning such as that illustrated in Figures 3-4. It would reduce the 
ram pressure of infall and the mass t o be photo-dissociated, makin g it easi e r to e j ect matter for a 
given core explosion mechanism. See lBrown. Bethe &: Bavml (|1982l ); iBethd (|1990l ); lArnettl (|1996l ). 



Changes in the mass and entropy of the collapsing co re will affec t its d ynamics. It has long 
been apparent that rotational effects too should play a role (|Hoyldll946l . ll955l ). The historical focus 
on non-rotating collapse was merely because the non-rotating problem was feasible with computing 
resources then available. All of these characteristics of the progenitor model might be modified 
significantly by the use of 3D initial models. The asymmetry and the rotational structure of the 
progenitor can be significantly modified by burning shell interactions, as can the ra te of f a llback , 
which depends upon the ma n tle structure and rotation a r ound the Fe core. See iFrverl (|1999); 
Fryer. Benz. &: Herantl (|1996l ); iFrver. Young fe Hungerfordl ( 20061 ) for a discussion of fallback and 
black hole formation. 

Core collapse is a converging flow (density increases); explosion is a diverging one (density 
decreases). A symmetries in convergent flow grow, which is a problem for inertial-confinement 
fusion efforts (|Lindllll998l ). In divergent slow, asymmetries decrease in importance, and the flow 
tends toward the spherically symmetric similarity solutions (|Sedovlll997l ). Because of the growth of 
asymmetry during collapse, it is important to have realistic estimates of the asymmetry in initial 
models of progenitors; ID mod els are spherical, so that the seeds of asy m metry are int r oduce d 
numerically or arbitrar ily (e.g., IWongwathanarat. Janka. fe Miillerl (|2010l ): IBruen. et al. I (|2009l ); 
Nordhaus. et al. I (|201ol) ). 
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The 2D simulations above suggest that progenitors prior to collapse develop large asymmetries 
in the O shell. What about the Fe core, which is what collapses? Si burning is active is a laye r 



of convective cells, so that the asymmetries will tend to average out (jArnett fc Meakinl l2011bl ). 
Simply scaling from Figure 0] suggests asymmetries in the Si shell of order of a few percent or 
more. In the absence of simulations including the Fe core in a dynamical way, we have no plausible 
quantitative estimates of asy mmetries in the Fe core i tself prior to collapse. Turbulence in O and Si 



shells will drive fluctuations (jArnett &i Meakinll2011bl ) ; the resulting motion will induce fluctuations 



i n the URCA shells in the core, and affect the change of entropy and electron fraction there (see 



(|Arnettlll996l ). §11 and §12). This in turn will affect the mass of the core at collapse, and thus the 
parameters of the collapse mechanism. Our understanding of the coupled physics of the dynamics 
of the core as it approaches collapse is still uncomfortably vague. 

Neutron Star Kicks. For core collapse progenitors, solitary or in binaries, the internal structure 
becomes increasingly inhomogeneous in radial density as e volution prec edes, and evolves toward a 



condensed core and extended mantle structure (Fig. 10.4 in lArnettl (119961 )). As the core plus mantle 
mass of these objects approaches the Chandrasekhar mass from above, this tendency increases. 
Supernovae of type SNIb and SNIc have light curves which require that they have such masses just 
prior to collapse. Consider a simple model in which there is a point-like core inside an extended 
mantle. If the core is not located at the center of mass, then the mantle must be displaced from the 
center of mass in the opposite direction. More of the mantle mass lies to one side of the core, so 
there is a net gravitational force which pulls the core back toward the center of mass (which does 
not move). Similarly the core exerts an equal and opposite pull on the mantle to bring the mantle 
back toward the center of mass. With no dissipation, an oscillation would ensue. The motion 
of the core relative to the fluid in the mantle generates waves in the mantle material, providing a 
means for dissipation, so that in the absence of driving, the oscillation of the core and mantle about 
the center of mass would be damped, and settle to a state in which both the core and mantle are 
centered on the center of mass. 

If there were a driving mechanism for core-mantle oscillation, there would be an asymmetry 
due to the displacement of core and mantle relative to the center of mass. The core collapse would 
give an off-center explosion within the mantle, even if the core collapse gave a perfectly spherical 
explosion shock relative to the center of mass of the core. 

Figures 3-4 above suggest that multiple shell burning may provide a driving mechanism for core- 
mantle oscillation; a computational domain containing an entire hemisphere would have allowed an 
£ = 1 mode to develop, which is even more suitable for driving such oscillations. Moreover, there 
will be a difference in the strength of driving depending upon the mass of the O burning shell; 
low mass shells will be less effective at driving the heavier core. Slow accretion onto an ONeMg 
white dwarf, or evolution to collapse by electron-capture, would occur by O burning in a low mass 
shell. However, evolution to core collapse by the instability of a more massive F e core generally 



occurs concurrently with more massive O burning shells, such as described above. Ivan den Heuvel 



(|2010l ) has suggested, on the basis of data on double neutron stars in the galaxy, that formation of 
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neutron stars from the collapse of ONeMg cores might occur with almost no kick velocity at birth, 
while neutro n stars formed by Fe core collapse wo uld receive a large space velocity at birth. See 
the review bv lKalogera. Valsecchi fe Willemsl (|2008l ) for background and references. The discussion 
above provides a physical mechanism for the empirical suggestion of van den Heuvel; the size of 
the kick velocity at collapse will depend upon the mass of the oxygen shell surrounding the core, 
and is driven by the dynamics of multiple shell burning. 



Wongwathanarat. Janka. fe Miillerl (|2010l ) find more vigorous kicks from collapse of Fe cores 
because the explosion develops sooner in ONeMg core, and the longer term instability in the post 
bounce behavior has less time to be effective. However, the collapse simulations to date have used 
small seed perturbations which may be unrealistic (they are much smaller than those found in 
Figures 3-4). 

In addition, the simulations shown in Figures 3-4 assume a static Fe core, and thus underesti- 
mate the total effect: asymmetry in the Fe core itself should have an important effect on the collapse 
and bounce, as mentioned above. A relatively massive mantle may itself affect the post-bounce 
behavior of the explosion by setting the initial condition whic h results in symmetry breaking (see 
below), which in the long term gives the hydrodynamical kick ([Wongwathanarat . Janka. fe Miiller 
2010 ). Either way, the higher kick velocity is associated with a more massive mantle, which provides 
the mechanism for symmetry breaking, and a complete picture needs to be developed. 

Early 7-rays. Gamma rays from the deca y of Ni 56 were observed in SN1987A before they were 
expected, based upon ID progenitor models (jArnett. Bahcall. Kirshner fe Wooslevll 19891 ). A strong 
O shell eruption, if hot enough to produce some Ni, followed by convective buoyancy and penetra- 
tive convection prior to collapse, would explain the early detection of gamma-rays, with no new 
hypotheses. Alternatively, explosive burning during the passage of the ejection shock would give 
an uneven distribution of fresh Ni 56 . Either way, some Ni 56 would be moved out further than in a 
spherically symmetric model, allowing earlier escape of 7-rays. 

Young SN Remnants. Young SN remnants have not yet mixed with the interstellar medium, 
and contain abundance information about the progenitor. The dynamic nature of pre-collapse evo- 
lution adds a new consideration to a ttempts to connect progenitor models to observations of young 
supernova remnants, s uch as Cas A (lYoung. et al. 1120081 ). For example, the puzzling inversion of Fe 



relative to Si found by lHughes. et al. I ([20001 ) in Cas A could easily be explained by vigorous dynam- 
ics of the O shell prior to core collapse. The spherical ID models are likely to be an inadequate ba- 
sis for interpreting observational data (e.g., Fesen. Becker fe Goodrich (jl988l ): Smith. J. P.. et al. 



(|2009l )). which may now be reanalyzed with a broader insight. Aspherical shock waves from the col- 
lapse d core become more spherical as they propagate outward ([Wongwathanarat. Janka. fe Miiller 
2010), so that even a very non-spherical collapse may be ineffective at producing asymmetries in 
the O shell. However, pre-existing asymmetries in the O shell, already significant, will be enhanced 
by explosive burning as the shock passes. 
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6. Summary 

The evolution of core-collapse progenitors is likely to be strongly dynamic, non-spherical, and 
may have extensive inter-shell mixing. These effects are ignored in existing progenitor models. 

The cyclonic patterns typical of 2D simulations are unstable in 3D, breaking apart and be- 
coming the turbulent cascade of Richardson and Kolmogorov. This enhances damping, and results 
in the lower velocities seen in 3D relative to 2D simulations. However, simulations in 3D will 
have essentially the sam e driving mechanisms as in 2D. Based upon existing simulations (e.g., 



Meakin fe Arnettl (|2007bl )). it is unlikely that the increased damping will eliminate dynamic behav- 
ior entirely. The increased damping may be able to moderate the eruptions seen in 2D, so that a 
set of quasi-steady dynamic pulses develops and continues until the core collapses. Alternatively, 
the increased damping may be inadequate to prevent continued growth of the instability, so that 
eruptions such as seen in Figure 2] will develop anyway, at a later time. The ultimate behavior 
would then be decided well into the non-linear regime. An extreme case would be an explosion 
powered by O burning prior to collapse of the Fe core. The observed light curve would depend 



upon the mass, kinetic energy and amount of ejected Ni 56 ( Arnettl 1 199a ). 



We need full 3D simulations to determine the quantitative impact of these new phenomena. 
From the discussion above it is possible to determine the features needed for such simulations: 

• Full 4ir steradians, including the whole core, to get the lowest order fluid modes [i = 1), 
rotation, and low order MHD modes, 

• Real EOS, to capture effect of electron-positron pairs and relativistic partial degeneracy, 

• Network, for realistic burning of C, Ne, O and especially Si with e-capture in a dynamic 
environment, 

• Multiple shells, (C, Ne, O, Si) to get shell interactions, 

• Sufficient resolution, to get turbulence and to calculate coherent structures (ILES), and 

• Compressible fluid dynamics, to get mixed mode waves, and possible eruptions. 



Low Mach number solvers such as MAESTRO (INonaka. et al. 1 120101 } may be useful, if gener- 
alized to include a dynamic background (the core evolution accelerates) , or applied to earlier stages 
of neutrino cooled evolution (which may be strongly subsonic). 

This is a challe nging combination of con straints, but such computations are becoming feasible. 



If we scale from the lMeakin k, Arnettl (|2007bl ) 3D simulation, the increased solid angle gives a factor 
of 85, the increase in radius a factor of 2, for a total increase of 170. A further increase in the radius 
of another factor of 2 would increase the computational load by a factor of 340, and would allow 
investigation of eruptions further into the strongly nonlinear regime than shown here. However that 
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simulation was done on a small Beowulf cluster (~ 100 cpus, which were slower than available now). 
This factor of 340 from the increased computational domain is more than balanced by the increased 
computational power available with top level machines. Doing an equivalent simulation, but in 3D, 
is feasible. More difficult is including the Fe core, which requires a dif ferent grid near the origin, but 
this h a s already been solved in different ways b y several groups (see [W oodw ard. Porter, fc Jacobs 
(|2003h : bearborn. Lattanzio. fc Egdetonl (|200fih : IWongwathanarat. Janka. fc Miillerl (|201ol ^. The 
computational demands, of a 3D simulation of the evolution of a progenitor into hydrodynamic 
core collapse, seems to be no worse than the computational demands of a single 3D core collapse 
calculation through bounce. 
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